% Readme
% 1.NOTE:Put question3.m and helicopter.m at the same path.
% 2.Structure: functions: helicopter.m
%              plotting: question3.m
%              answers: question3.m 
% 3.Assumption: 
%    linear twist; 
%    Uniform inflow velocity;
%    no taper; the induced velocity is constant along the blade
%    the chord length is constant
%    The Airfoil lift and drag characteristics are the same as for NACA0012
%    The rotor is far from the ground.
% 4.User guide
%  4.1 Change the parameters 
%
%    Change the basic parameters, go to properties parts of class
%      helicopter.m and retype them.
%    Change plotting parameters, go to question3.m Question i-viii retype
%      them.
%    All the methods of question i-viii are shown below, the parameters can be easily
%    changed
%    
%  4.2 Add the parameters
%    Go to properties parts of class helicopter.m, just add a parameter.
%    When using it, using helicopter.parameter. 
clear;clc;

N=1000;
mu_s        = zeros(1,N+1); % matrix of mu
power_r     = zeros(1,N+1); % matrix of required power
power_a     = zeros(1,N+1); % matrix of available
power_e     = zeros(1,N+1); % matrix of excess power
v_f         = zeros(1,N+1); % matrix of forward speed
V_c         = zeros(1,N+1); % matrix of climb speed
Cp          = zeros(1,N+1); % matrix of Cp
theta_climb = zeros(1,N+1); % matrix of climb angle

mu_max      = 0.45; % maximum of mu
delta_mu    = mu_max/N; %step of mu
mu          = 1e-8; % to avoid getting NAN

for i=1:N+1
    mu_s(1,i) = mu;
    power_r(1,i) = helicopter.power_required(mu); 
    power_a(1,i) = 1544;    % individual parameters 
    power_e(1,i) = power_a(1,i) - power_r(1,i);
    V_c(1,i) = helicopter.climbRate(mu);
    theta_climb(1,i) = helicopter.theta_climb(mu);
    v_f(1,i) = helicopter.omega * helicopter.R * mu;
    Cp(1,i) = helicopter.Cp(mu);
    mu = mu + delta_mu; % start next iteration
end

% Question i
% Plot the figure of power required and power available
figure(1);
plot(mu_s, power_r, 'r', 'LineWidth', 2); hold on;
plot(mu_s, power_a, 'b', 'LineWidth', 2); 
legend('power required','power available', 'Location','SouthWest');
xlim([0,mu_max]);
ylim([0, 2000]);
xlabel('\mu','FontSize',11);
ylabel('power/kw', 'FontSize', 11);
title('power required and available at sea level vs. advance ratio');

% question ii
% plot the node of power required and power available
figure(2);
plot(mu_s, power_r, 'r', 'LineWidth', 1); hold on;grid on; grid minor; % grid
plot(mu_s, power_a, 'b', 'LineWidth', 1); 
k = find(abs(power_r-power_a)<=0.5); % using find function to get the node
plot(mu_s(k),power_r(k),'*');
text(mu_s(k)-0.1,power_r(k),['(',num2str(mu_s(k)),',',num2str(power_r(k)),')']); % to show text inside the graph
legend('power required','power available','max speed', 'Location','SouthWest');
xlim([0,mu_max]);
ylim([0, 2000]);
xlabel('\mu','FontSize',11);
ylabel('power/kw', 'FontSize', 11);hold off;
title('The maximum forward speed at sea level');
v_max = helicopter.omega * helicopter.R * mu_s(k) % maximum of forward flight

% question iii
% plot the figure of power excess
figure(3);
plot(mu_s,power_e,'r','LineWidth', 2);
% xlim([0,mu_max]);
ylim([0, 1050]);
xlabel('\mu','FontSize',11);
ylabel('power/kw', 'FontSize', 11);
legend('power excess', 'Location','SouthWest');
title('excess power at sea level vs. advance ratio');

% question iv
% plot the climb rate of the helicopter
figure(4);
[a1,a2] = max(V_c); % a1 output the maximum value, a2 returns its order number.
max_climb = V_c(a2);
plot(mu_s(1,:),V_c,'r','LineWidth', 2);hold on;
plot(mu_s(a2),max_climb,'*');
xlim([0,mu_max]);
% ylim([0, 2000]);
xlabel('\mu','FontSize',11);
ylabel('rate of climb m/s', 'FontSize', 11);hold off;
legend('power excess', 'max climb rate','Location','SouthWest');
text(mu_s(a2),V_c(a2),['(',num2str(mu_s(a2)),',',num2str(V_c(a2)),')']);
title(' rate of climb at sea level vs. advance ratio');

% question v
% calculate the climb angle
figure(5)
plot(mu_s(1,2:N+1),theta_climb(1,2:N+1),'r','LineWidth', 2);
legend('climb angle','Location','SouthWest');
ylim([0, 1.6]);
title('climb angle vs advance ratio');
xlabel('\mu','FontSize',11);
ylabel('angel/rad', 'FontSize', 11);

% question vi
% calculate the minimum power speed
figure(6);
plot(mu_s, power_r, 'r', 'LineWidth', 1);hold on;grid on; grid minor;
[b1,b2] = min(power_r); % b1 output the minimum value, b2 returns its order number.
plot(mu_s(b2),power_r(b2),'*');
text(mu_s(b2),power_r(b2),['(',num2str(mu_s(b2)),',',num2str(power_r(b2)),')']);
xlabel('\mu','FontSize',11);
ylabel('power/kw', 'FontSize', 11);
title('power required and available at sea level vs. advance ratio');
v_min = helicopter.omega * helicopter.R * mu_s(b2)

% question vii
% calculate the endurance
v_maxEndurance = mu_s(b2) * helicopter.omega * helicopter.R;
WF = (helicopter.m_empty + 0.1 * helicopter.m_maxFuel + helicopter.m_payLoad)*9.81;
W0 = (helicopter.m_empty + helicopter.m_fuelLoad + helicopter.m_payLoad)*9.81; %initial helicopter weight
E_endurance = helicopter.thrust(b2)/(helicopter.consu*1544*1000)*log(W0/WF);

% question viii
% calculate the range
figure(8);
plot(mu_s,power_r, 'r', 'LineWidth', 2);hold on;grid on; grid minor;
c1 = find(abs(Cp/mu_s - gradient(Cp)./gradient(mu_s))<=0.000004);
% using find and gradient function to get the mu that satisfy Cp/mu=dCp/dmu
plot(mu_s(c1),power_r(c1),'*');
text(mu_s(c1),power_r(c1),['(',num2str(mu_s(c1)),',',num2str(power_r(c1)),')']);
xlabel('\mu','FontSize',11);
ylabel('power/kw', 'FontSize', 11);
title('power required and available at sea level vs. advance ratio');

v_maxRange = mu_s(c1) * helicopter.omega * helicopter.R;
R_max = helicopter.thrust(c1)*v_maxRange/(helicopter.consu*1544*1000)*log(W0/WF)
